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ABSTRACT 

We use 5000 cosmological A-body simulations of lft,~ 3 Gpc 3 box for the concordance ACDM model 
in order to study the sampling variances of nonlinear matter power spectrum. We show that the 
non-Gaussian errors can be important even on large length scales relevant for baryon acoustic oscil- 
lations (BAO). Our findings are (1) the non-Gaussian errors degrade the cumulative signal-to-noise 
ratios (S/N) for the power spectrum amplitude by up to a factor of 2 and 4 for redshifts z = 1 
and 0, respectively. (2) There is little information on the power spectrum amplitudes in the quasi- 
nonlinear regime, confirming the previous results. (3) The distribution of power spectrum estimators 
at BAO scales, among the realizations, is well approximated by a Gaussian distribution with variance 
that is given by the diagonal covariance component. (4) For the redshift-space power spectrum, the 
degradation in S/N by non-Gaussian errors is mitigated due to nonlinear redshift distortions. (5) 
For an actual galaxy survey, the additional shot noise contamination compromises the cosmological 
information inherent in the galaxy power spectrum, but also mitigates the impact of non-Gaussian 
errors. The S/N is degraded by up to 30% for a WFMOS-type survey. (6) The finite survey volume 
causes additional non-Gaussian errors via the correlations of long-wavelength fluctuations with the 
fluctuations we want to measure, further degrading the S/N values by about 30% even at high redshift 
z = 3. 

Subject headings: cosmology: theory - large-scale structure of universe 



1. INTRODUCTION 

Baryon acoustic oscillations (BAO) are imprinted in 
the distribution of galaxies, of which the characteris- 
tic length scale can be used as a standard ruler in the 
universe (e.g. Eisenstein, Hu & Tegmark 1998; Blake 
& Glazebrook 2003; Seo & Eisenstein 2003; Matsubara 
2004; Guzik, Bernstein k, Smith 2007). BAO provides a 
powerful way of probing the nature of dark energy. Large 
galaxy surveys such as the Sloan Digital Sky Survey and 
two degree Field Survey detected the BAO signature and 
provided constraints on the dark energy equation of state 
(Cole et al. 2005; Eisenstein et al. 2005; Percival et al. 
2007; Okumura et al. 2008; Gaztanaga, Cabre & Hui 
2008; Sanchez et al. 2009). Future larger surveys are 
aimed at measuring the BAO scale more accurately and 
hence yielding tighter constraints on the nature of dark 
energy (see e.g. Benitez et al. 2008). 

The BAO signature in the galaxy power spectrum is 
very small, of the order of a few percent modulation 
in amplitude, and hence measurements of the precise 
length scale are hampered by a number of effects. For ex- 
ample, nonlinear gravitational evolution, redshift space 
distortion, galaxy formation processes and the associ- 
ated scale-dependent bias, all compromise a robust de- 
tection. Accurate theoretical models are clearly needed. 
A number of authors resort to using numerical simula- 
tions (Meiksin, White & Peacock 1999; Seo & Eisenstein 
2005; Huff et al. 2007; Smith, Scoccimarro & Sheth 
2007, 2008; Angulo et al. 2008; Takahashi et al. 2008; 
Seo et al. 2008; Nishimichi et al. 2008) whereas oth- 



ers use perturbation theory (Crocce & Scoccimarro 2006, 
2008; Jeong & Komatsu 2006, 2009; Nishimichi et al. 
2007; McDonald 2007; Matarrese & Pietroni 2007, 2008; 
Pietroni 2008; Matsubara 2008a, b; Taruya & Hiramatsu 
2008; Takahashi 2008; Nomura, Yamamoto & Nishimichi 
2008; Rassat et al. 2008). It is important to note that 
one needs accurate estimates not only for the power spec- 
trum but also its covariance (e.g. Scoccimarro, Zaldar- 
riaga & Hui 1999; Meiksin & White 1999; Habib et al. 
2007). The covariance describes statistical uncertainties 
of the power spectrum measurement as well as the band 
powers at different wavenumbers are correlated with each 
other. Hence once the well-calibrated covariance is ob- 
tained, one can derive unbiased, robust constraints on 
cosmological parameter from the measured power spec- 
trum (see Ichiki et al. 2008 for such an example to show 
the importance of the covariance estimation). 

The power spectrum covariance matrix has only di- 
agonal elements for the Gaussian density fluctuations 
(e.g. Feldman, Kaiser & Peacock 1994). The relative 
error of the power spectrum of a given wavenumber is 
then simply given by the square root of the number of 
Fourier modes available from the survey volume. How- 
ever, at small length scales, non-vanishing off-diagonal 
parts of the covariance arise due to the mode coupling 
(Scoccimarro, Zaldarriaga & Hui 1999; Meiksin & White 
1999; Smith 2008). This non-Gaussian contribution is 
described by the trispectrum or the Fourier transform of 
the 4-point correlation function. Cooray & Hu (2001) 
used the halo model to estimate the trispectrum contri- 
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bution and showed that the non-Gaussian errors do de- 
grade the precision of cosmological parameter determina- 
tion, and therefore cannot be ignored for planned future 
surveys (see also Takada & Jain 2008; Eifler, Schneider 
& Hartlap 2008). Also recently, Smith (2008) studied the 
covariance matrix of the halo power spectrum using nu- 
merical simulations. Sefusatti et al. (2006) also studied 
the power spectrum covariance using PTHalos (Scocci- 
marro & Sheth 2002). 

In this paper, we use an unprecedentedly large num- 
ber of simulation realizations to estimate the covariance 
matrix of the matter power spectrum in both real and 
redshift space. Our sample is more than 2 orders of 
magnitude larger than those used in the previous works, 
yielding well-converged estimates on the power spectrum 
covariance. We compare our simulation results with the 
analytical estimates based on perturbation theory and 
halo model. In these comparisons we also include the 
new effect of non-Gaussian errors that inevitably arise 
for a finite-volume survey, as first pointed out in Rimes 
& Hamilton (2005; also see Rimes & Hamilton 2006; 
Hamilton, Rimes & Scoccimarro 2006; Neyrinck, Sza- 
pudi & Rimes 2006; Neyrinck & Szapudi 2007, and Lee 
& Pen 2008 for the observational implication based on 
the SDSS data). By using this large number of the real- 
izations, we also study how the power spectrum estimates 
are distributed in different realizations, i.e. the probabil- 
ity distribution of power spectrum, and then compute 
the higher-order moments, skewness and kurtosis, to ex- 
amine the overall impact of power spectra at high sigma 
ends. Furthermore we estimate the expected signal-to- 
noise ratio for measuring the power spectrum for a future 
galaxy survey, taking into account the shot noise contam- 
ination and the non-Gaussian errors. 

Throughout the present paper, we adopt the standard 
ACDM model with matter density £l m = 0.238, baryon 
density f2b = 0.041, cosmological constant tt\ = 0.762, 
spectral index n s — 0.958, amplitude of fluctuations 
as = 0.76, and expansion rate at the present time 
Hq = 73.2km s _1 Mpc" 1 , which are consistent with the 
WMAP 3-year results (Spergel et al. 2007). 

2. NUMERICAL SIMULATIONS 

We use the cosmological simulation code Gadget-2 
(Springel, Yoshida & White 2001; Springel 2005). We 
employ 256 3 particles in a volume of lOOO/i -1 Mpc on a 
side. We generate initial conditions of the seed density 
perturbations at z = 20 based on the standard Zel'dovich 
approximation using the matter transfer function calcu- 
lated by CAMB (Code for Anisotropies in the Microwave 
Background: Lewis, Challinor & Lasenby 2000). We ran 
5000 realizations of Particle Mesh (PM) simulations for 
the fiducial cosmological model, and use the snapshot 
outputs at z — 3,1 and to study the power spectrum 
covariances. 

To calculate the Fourier transform of the density field, 
denoted as <5(k), we first assign the A-body particles onto 
A| rid = 512 3 grids based on the cloud-in-cell method 
and then perform FFT 1 . We also correct the effect of 
the cloud-in-cell assignment scheme as S(k) — » 8(k) x 
[sinc(fc x i/2A gri d) sinc(fc y L/2A gri d) sinc(fc z L/2A gr i d )]~ 2 

1 FFTW home page: http://www.fftw.org/ 



with sinc(x) = sinx/x (Hockney & Eastwood 1988; An- 
gulo et al. 2008). The binned power spectrum for a given 
realization is estimated as 



P(k) 
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where the summation runs over all the Fourier modes 
whose length is in the range k — Ak/2 < |k| < k + Ak/2 
for a given bin width A A;. Here A& is the number of 
modes taken for the summation and is given as N k — 
E|k|efc ~ 47rfc 2 Afc/(2 7 r/ J L) 3 = Ffc 2 Afc/(27r 2 ) for the 
limit k 3> 1/i, where L is the simulation box size and 
V is the volume given by V — L 3 . The shot noise is 
not subtracted, since this effect is very small. The en- 
semble average of the power spectrum estimator is then 
computed by averaging the estimated spectra over the 
realizations: P(k) = {P(k)}. 

We have checked that our simulation result for 
the power spectrum agrees with the higher resolution 
TreePM result within 1%(3%) at k < 0.2(0.4)/i/Mpc 2 
(here the Nyquist wavenumber is k = 0.8/i/Mpc). If the 
initial redshift is set to be higher, e.g. z = 50, the re- 
sults agree within 2% for k < 0.2ft,/Mpc and 10% for 
k < 0.4ft,/Mpc for z = 0, 1,3. This is sufficient for our 
purpose, which is to estimate the impact of nonlinear 
clustering on the power spectrum covariances at BAO 
scales. 

3. COVARIANCE MATRIX 

The covariance between the power spectra, P{k\) and 
P{k2), is estimated from the simulation realizations and 
can be formally expressed in terms of the Gaussian and 
non-Gaussian contributions (e.g. Scoccimarro, Zaldar- 
riaga & Hui 1999; Meiksin & White 1999): 



cov(fc 1 , fc 2 ) EE ( P(fci) - P(ki) P(k 2 ) - P{k 2 ) 



+ 



V 



d 3 k[ d 3 k' 2 

met, J\W 2 \ek 2 Vki Vk 2 



T(k' 1 ,-ki,k^-k^),(2) 



where T is the trispectrum, and 8j? ik2 is the Kronecker- 
type delta function defined such that 8f? i1t , 2 = 1 if ki =k 2 
within the bin width, otherwise zero. The integration 
range in the second term is, as in Eq. (1), confined to 
the Fourier modes lying in the range k\ — Ak/2 < k < 
ki + Ak/2, and V? £i (i — 1,2) denotes the integration 
volume in Fourier space given by Vu 1 ~ 47rfc 2 Afc for the 
case of k 3> Afc. 

The first term of the covariance matrix represents the 
Gaussian error contribution ensuring that the two power 
spectra of different wavenumbers are uncorrelated, while 
the second term gives the non-Gaussian errors that in- 
clude correlations between power spectra at different fc's 
arising from nonlinear mode coupling. Both the terms 
scale with the simulation box volume as oc 1 /V . It should 
be also noted that the non-Gaussian term does not de- 
pend on the bin width (because J| k /i efc rf 3 k'/Vfc 1), 

2 The agreement is achieved in real space. In redshift space. 
PM simulations somewhat underestimate the power spectrum by 
20(10)% at z = 0, 1(3) at small length scales (k = 0.4h/Mpc). 
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so increasing Afc only reduces the Gaussian contribution 
via the dependence Nk oc Afc. However, the cumulative 
signal-to-noise ratio we will study below is independent 
of the assumed Afc. 

We will compare the simulation results with two an- 
alytical approaches to estimate the covariance matrix: 
(1) perturbation theory and (2) halo model. In per- 
turbation theory, following Scoccimarro et al. (1999; 
also see Neyrinck & Szapudi 2008), the power spectrum 
and trispectrum are, self-consistently including up to the 
third order perturbations of <5, expressed as 

P 2 {k 1 ) = Fg,{kx) + 2flin(fci) [P 22 (fci) + Fi 3 (fci)] , 
T(k 1 ,-k 1 ,k 2) -k 2 ) = 12 Pr m (fci)Pi„(fc 2 ) 
x [ J F3(k 1) -k 1 ,k 2 )fl in (fc 1 ) + (fci «-» k 2 )} 
+8P Iin (|ki - k 2 |) [P 2 (k x - k 2 ,k 2 )P lin (fc 2 ) 



+ (fci «-»fc)f , 



(3) 



where Pn n denotes the linear-order spectrum, and P 22 
and Pi 3 are the one-loop corrections to the nonlinear 
power spectrum (Makino, Sasaki & Suto 1992; Jain & 
Bertschinger 1994) and P 2 and P3 are the kernels for 
the second and third order density perturbations (e.g. 
Bernardeau et al. 2002). 

In the halo model, the power spectrum is given by a 
sum of two terms, the so-called one-halo term and two- 
halo term (Seljak 2000; Ma & Fry 2000; Peacock & Smith 
2000; also see Cooray & Sheth 2002 for a review). Simi- 
larly, the trispectrum consists of four terms, from one to 
four halo terms: 



T = T 



ih 



T 



2h 



T 



T 



■4h 



(4) 



The explicit expressions of each term can be found in 
Cooray & Hu (2001). In nonlinear regime, i.e. large k, 
the 1-halo term gives dominant contribution to the total 
power of trispectrum, while the different halos terms be- 
come more significant with decreasing k and the 4-halo 
term that includes the PT trispectrum contribution be- 
comes dominant on very small k. 

To complete the halo model approach, we need suitable 
models for the three ingredients: the halo mass function 
(Shcth & Tormen 1999), the halo bias parameters (Mo & 
White 1996; Mo, Jing & White 1997), and the halo mass 
density profile (Navarro, Frenk & White 1997), each of 
which is specified by halo mass m and redshift z for a 
given cosmological model. The details of our halo model 
implementation can be found in Takada & Jain (2003; 
2008). 

In our previous paper (Takahashi et al. 2008), we 
found that a finite size simulation causes the growth of 
large-scale density perturbations to be deviated from the 
linear theory prediction, and the deviation is well de- 
scribed by the nonlinear mode coupling. We proposed a 
method to "correct" the deviation in the finite size sim- 
ulations to obtain the ensemble averaged expectation of 
power spectrum without running ideal simulations with 
infinite volume (practically very large volume), as also 
demonstrated in Nishimichi et al. (2008). In this paper, 
where we are discussing the power spectrum covariance 
for some finite survey volume, we do not have to "cor- 
rect" the deviation in power spectrum of each realization, 
because the scatters are already included in the covari- 
ance formula (the terms with P 2 in Eq. [3]). 



We use 5000 realizations of each output redshift to 
directly estimate the covariance matrix according to 
Eq. @. To be more explicit, denoting the power spec- 
trum of the i-th realization as Pj(fc), we can estimate the 
covariance as 



cov(fci, fc 2 ) 



1 



r i=l 



[Pi(h)-P(h)}[PM-P(k 2 )}, 

(5) 

where N T is the number of realizations, i.e. N T — 5000 in 
our case, and P(fc) denotes the mean spectrum computed 
as P(k) = (1/iVr) Pi(k)- As shown in Appendix, the 
accuracy in estimating the covariances scales with the 
number of realizations used 3 . For example, the relative 
accuracy of estimating the diagonal covariance elements 
is found to scale approximately as (iV r /2) _1 ' 2 . Hence, 
with the aid of 5000 realizations, we can achieve a few 
%- level accuracies in estimating each elements of the co- 
variance, an improvement by an order of magnitude over 
previous works. 

4. COMPARISON WITH THEORETICAL MODELS 
4.1. Results in Real Space 

Fig[T] shows the diagonal elements of the covariance 
matrix as a function of wavenumbers. The diagonal el- 
ements plotted are divided by the Gaussian covariances 
of linear power spectra at each redshift such that the 
values become unity in the linear regime limit (k — ► 0). 
Therefore the deviations from unity arise from the non- 
linear evolution of P(fc) and the non-Gaussian covari- 
ance contribution. The cross, triangle and circle sym- 
bols show the simulation results for redshifts z = 3, 1 
and 0, respectively. Note that we adopt the bin width 
of Afc = 0.01/iMpc -1 throughout this paper. The de- 
viations from the Gaussian errors become more signifi- 
cant at lower redshifts. For comparison, the solid curves 
show the analytical predictions obtained when the per- 
turbation theory is employed to estimate the covariances 
as described around Eq. The perturbation theory 
fairly well reproduces the simulation results within 20% 
up to k < 0.24/i/Mpc at z = and k < 0.4/i/Mpc at 
z = 1 and 3, respectively. However, at lowest redshift 
z — 0, stronger nonlinear effects are seen even on these 
large length scales corresponding to the BAO scales. The 
dashed curves show the halo model results which take 
into account of this nonlinear effect. The halo model 
fairly well fits the simulation results over the range of 
wavenumbers studied. At z — 1,3 the PT predicts the 
larger variance than the halo model, because the one-loop 
power spectrum in the Gaussian term ^ overestimate 
the power spectrum. 

FigE] shows the off-diagonal elements of the covariance 
matrix. For illustrative purpose we study the correlation 
coefficient matrix defined as 



r(fci, fc 2 ) = 



cov(fci, fc 2 



(6) 



i/cov(fci, fci)cov(fc 2 , fc 2 ) 

The coefficients are normalized so that r = 1 for the 
diagonal components with k\ = fc 2 . For the off-diagonal 

3 More precisely, the accuracy of estimating the covariance is 
determined by the covariance of the power spectrum covariance 
that includes up to the 8-point correlation functions (Kayo et al. 
in preparation). 
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Fig. 1. — The diagonal components of the power spectrum co- 
variance as a function of wavenumbers, for redshifts z = 0, 1 and 3. 
The results are divided by the Gaussian covariance of the linearly 
evolving power spectrum. Therefore the deviations from unity arise 
from both the nonlinear clustering and the non-Gaussian errors. 
The symbols are the simulation results, while the solid curves show 
the results obtained when the perturbation theory (PT) is used to 
compute the non-Gaussian covariance. The dashed curves show 
the halo model results. 

components r — > 1 implies strong correlation between the 
two spectra, while r = corresponds to no correlation. 
Note again that the matrix elements r depend on the 
bin width: a finer binning, i.e. a smaller Ak 7 decreases 
the off-diagonal components. First of all, comparing the 
three panels of Fig.[2]mamfests that the off-diagonal com- 
ponents have greater amplitudes with increasing k. The 
PT fits the data within 0.1 for k < 0.15/i/Mpc at z = 0, 
for k < 0.24/i/Mpc at z = 1 and for k < 0.29/i/Mpc 
at z = 3. For the redshift dependence, there is almost 
no cross-correlations at redshift z = 3, while there are 
increasing cross-correlations at lower redshifts. 

The PT results start to underestimate the correlation 
strengths with increasing k and at lower redshifts due to 
the stronger nonlinearities. Compared to Fig. [TJ the PT 
results are found to be less accurate to describe the off- 
diagonal components at z = 1 and 0. The dashed curves 
are the halo model results, which are in a good agreement 
with the simulation results, especially at z = 0. We 
found that an inclusion of the 2- and 3-halo terms is 
important to describe the scale dependences of the off- 
diagonal correlations. However, the halo model displays 
a sizable disagreement at some scales, and is not well 
accurate. Therefore a further refinement of the model 
predictions based on this kinds of large-scale simulations 
is needed to accurately model the measurement errors 
of power spectrum, especially for future high-precision 
surveys. 

4.2. Results in Redshift Space 

In this section, we examine the covariance of the 
redshift-space power spectrum that is a more direct ob- 
servable in galaxy surveys. The redshift-space power 
spectrum in each realization is computed as follows. As- 
suming the distant observer approximation, we first cal- 
culate the density perturbations in redshift space as de- 
scribed in § [2] but properly taking into account modula- 



tions of iV-body particle positions in redshift space due 
to the peculiar velocities. The density perturbation field 
is thus given as a function of wavenumbers k\\ and 
that are parallel and perpendicular to the line-of-sight 
(taken from one direction in the simulation box). As a 
result, the redshift-space power spectrum P s is given as a 
two-dimensional function due to the statistical isotropy: 
P s (kj_,k\\). In this paper, for simplicity, we focus on 
the spherically averaged redshift-space spectrum over the 
shell of a radius k with the width Afc: 



* |k'|efe 



lMMh k l)l 5 



(7) 



where k' = y fcj 2 + |k^_| 2 and Nk is the number of modes 

in the spherical shell in redshift space. Likewise, the co- 
variance matrix of P s n can be estimated by averaging the 
spectrum estimators among the simulation realizations 
as in Eq. ©. 

According to the linear perturbation theory of struc- 
ture formation, the redshift-space power spectrum can 
be simply related to the real-space spectrum under the 
distant observer approximation as P s {ku,k±) — (1 + 
//z 2 ) 2 P(/c), where \x = k\\/k is the cosine between the 
line-of-sight and the wavevector and / is the linear red- 
shift distortion, expressed in terms of the linear growth 
rate Di as / = {d\xiDx/d\xia)/b (Kaiser 1987) with the 
bias parameter b = 1 for the dark matter power spec- 
trum. Note that all the spectra we have considered are 
for the total matter distribution. Averaging the redshift- 
space spectrum over the cosine angle fi yields the linear 
theory prediction that is to be compared with the simu- 
lation result given by Eq. ([7]): 



P s0 (fc) = [l + (2/3)/ + (l/5)/ 2 ] P(k). 



(8) 



The prefactor in front of P(k) on the r.h.s. of Eq. ((8) 
does not depend on wavevector. Hence, from Eqs. Ij2|) 
and ijSj), the linear theory tells that the covariance of 
the redshift-space power spectrum (jHJ) can be simply ex- 
pressed as 4 

cov s (A;i,fc 2 ) = /fp s o(fci)-P s o(fci)') (Ps Q {k 2 )- P sQ {k 2 ) 



|/ 3 + |/ 4 



1 y ki 



§/+i/ 2 



i+lz+l/ 2 )' 



(9) 



Due to the additional factor that depends solely on /, the 
covariance amplitude of P s o is greater than the standard 
Gaussian error, (2/iV fe )P 2 , by 6, 16 and 20% at z = 0, 1 
and 3 for the ACDM model, respectively. Note that the 
covariance form (J5|) is valid only for the asymptotic limit 
of large length scales, and in general the nonlinear clus- 
tering effects cause deviations from the Kaiser formula 
on the BAO scales (e.g., Scoccimarro 2004). 

The diamond symbols in FigJ3]show the diagonal com- 
ponents of the redshift-space power spectrum covariances 
as a function of wavenumbers and at three output red- 
shifts. The diagonal components are divided by the stan- 
dard Gaussian errors, 2/./VfcP 2 (fc), where we have used 

4 The angular average of the covariance is proportional to 
f dfj,(l + //Lt ) , while the square of the angular averaged P s (k) 
is proportional to [ f dfi(l + ffj, 2 ) 2 ] 2 . Hence the two quantities are 
not same and the extra factor in Eq. ([9jl appears. 



Covariance matrix of the matter power spectrum 



5 



k 1= 0.035h/Mpc 



0- *» 



real space 
Ak=0. 01 h/Mpc 



z=3 
z=1 
z=0 




- 


k 1= 0.135h/Mpc 




PT 




HM 











k,=0.235h/Mpc 




0.2 



0.4 



0.2 



0.4 



0.2 



0.4 



k 2 (h/Mpc) 



k 2 (h/Mpc) 



k 2 (h/Mpc) 



Fig. 2. — The correlation coefficient matrix r(ki,k2), defined in Eq. JB), as function of k%, where fci is kept fixed to k\ = 0.035, 0.f35 
and 0.235/i/Mpc in the left, middle and right panels, respectively. The solid curves denote the PT predictions, while the dashed curves 
show the halo model results. The simulation results at z = and 1 show greater amplitudes in the off-diagonal covariances than the PT 
predictions. 
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Fig. 3. — The diagonal covariance components as a function of wavenumbers, for real- and redshift-space power spectra. Note that the 
redshift-space power spectrum studied here is the spherically averaged spectrum over a shell of a given radius k in redshift space (see text 
for the details). We show the covariances divided by the Gaussian error contribution (the first term in Eq. [2]): at large length scale limit 
[k — * 0), the real- and redshift-space values approach to unity (solid line) and to the constant factor that is given by the Kaiser's linear 
distortion (dashed line), respectively. The non-Gaussian error contribution is relatively suppressed in redshift space due to the nonlinear 
redshift distortions. 




ki (h/Mpc) ki (h/Mpc) 

Fig. 4. — The correlation matrix at z = (left panel) and 2 = 1 (right panel). In each panel, the upper-left matrix elements are the 
off-diagonal covariances in real space, while the lower-right elements are for redshift space. 
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Fig. 5. — Probability distribution of the power spectrum es- 
timators P in the 5000 realizations. The solid curves show the 
Gaussian distribution with zero mean and variance that is set to 
the diagonal covariance component measured from the simulations 
at a given wavelength. 

the nonlinear power spectrum measured from the simu- 
lations. Note the difference in the normalization factor 
from that in FigfT] The horizontal line in each panel 
shows the prefactor in Eq. ([9]), the amplification fac- 
tor expected from the Kaiser's formula at the large-scale 
limit. Therefore the deviations from the horizontal line 
may come from two contributions: (1) the non-Gaussian 
error contribution caused by nonlinear clustering and 
(2) the nonlinear redshift distortions such as the effect 
caused by the virial motions within and among halos, 
known as the finger-of-God effect. The circle symbols 
denote the simulation results for the real-space spectrum 
computed in a consistent way, i.e. divided by the non- 
linear spectrum. The nonlinear effects on the covariance 
become more significant with increasing wavenumber and 
at lower redshifts. Interestingly, however, comparing the 
real- and redshift-space results manifests that the relative 
importance of the non-Gaussian covariances is weaker in 
redshift space, implying that the finger-of-God redshift 
distortions at small length scales more preferentially sup- 
press the covariance amplitudes than the power spectrum 
amplitudes (also see Meiksin & White 1999). 

Figf4] shows both the off-diagonal components of the 
covariances in real space (the left-upper elements in each 
panel) and redshift space (right-lower) , at redshifts z = 
and z — 1. The cross-correlations are more signifi- 
cant with increasing wavenumbers, while the correlation 
strengths are relatively weaker in redshift space. 

5. PROBABILITY DISTRIBUTION OF THE POWER 
SPECTRUM ESTIMATOR 

We have so far discussed the non-Gaussian covariance 
of the power spectrum estimator P. It would be also 
intriguing to study how the nonlinear clustering causes 
a non-Gaussian distribution in the power spectrum es- 
timators of a given k among our 5000 realizations. For 
example, if the estimators have a skewed distribution, a 
prior knowledge on the full distribution may be needed 
to obtain an unbiased estimate on the ensemble averaged 
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Fig. 6. — The skewness (top panel) and the kurtosis (bottom 
panel) of the power spectrum distribution shown in the previous 
figure, as a function of k for outputs at z = 20 and 0. The cir- 
cles and crosses are the simulation results, while the solid curves 
are the theoretical predictions expected when the power spectrum 
X 2 



band power at each k from a small number of realizations 
or a finite volume survey. Note that the power spectrum 
covariance simply reflects the width (variance) of the full 
distribution at each k's but does not contain full infor- 
mation on the probability distribution. 

Fig. [5] shows the probability distribution of the power 
spectrum estimators P among 5000 realizations, where 
we mean by "probability" that the distribution is normal- 
ized so as to give unity if the distribution is integrated 
over the x-axis values (see below). The cross, triangle, 
and circle symbols show the results for k = 0.05, 0.2, and 
0.4/i/Mpc, respectively. The distribution is plotted as a 
function of (N k /2) 1 / 2 (P / P - 1) for each k such that the 
mean and variance of the distribution are equal to zero 
and unity when the power spectrum distribution obeys 
the linear-regime Gaussian distribution. The simulation 
results show that the distribution is broadened with in- 
creasing k due to the stronger nonlinearities. The solid 
curves show the expected Gaussian distribution where its 
variance is set to the diagonal covariance measured from 
the simulations at each k, i.e. the variance includes the 
non-Gaussian covariance contribution as given in Fig. [3l 
Interestingly, the simulation results are rather well ap- 
proximated by the Gaussian distribution even in the non- 
linear regime. 

The remaining small deviations from the Gaussian dis- 
tribution can be quantified by studying the skewness 5*3 
and kurtosis S4 defined as 



S 3 = 



((P(fc)-P(fc)) 3 ) 
((P(k) -P(fc)) 2 ) 3 / 2 ' 
= ((P(fc)-F(fc)) 4 ) 
4 ((P(fe)-P(fe)) 2 ) 2 



(10) 
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The 53 and S4 are vanishing for the Gaussian distri- 
bution. If the density field obeys the random Gaus- 
sian fields, which is a good approximation in the linear 
regime, the power spectrum estimator of a given k (see 
Eq. Q]) obeys the XN k -distribution in analogy with the 
CMB power spectrum (Knox 1995). In this case, as de- 
rived in Appendix B, the skewness and kurtosis can be 
analytically computed as 

Fig. [5] shows the simulation results for S3 and S4 as a 
function of k at z = 20 and 0. The results are for a vol- 
ume of V = l/i~ 3 Gpc 3 , and the Sa t 4, scale as S3 oc y -1 / 2 
and S4 oc V~ x from Eq.fTTjl. The solid curves are the 
theoretical predictions of Eq. lfTTj) which well match the 
simulation results. Note that the skewness is positive, be- 
cause the x 2 -distribution has a long tail at large ends of 
P. Both 5*3 and S4 asymptotes to zero at high k, i.e., the 
probability distribution approaches to a Gaussian distri- 
bution at high k. The skewness grows from z — 20 to 
through the non-linear gravitational evolution, however, 
its value (£3 < 0.1) is very small. 

6. EFFECTS OF NON-GAUSSIAN COVARIANCE ON 
SIGNAL-TO-NOISE RATIO 

A useful way to quantify the impact of the non- 
Gaussian errors is to study the cumulative signal-to-noise 
ratio (S/N) for measuring the power spectrum over a 
range of wavenumbers, which is also sometimes called the 
Fisher information content (e.g. Tegmark et al. 1997). 
The S/N is defined, using the covariance as 

(|) 2 = S ^cov- 1 ^^)/^), (12) 

where cov -1 denotes the inverse of the covariance matrix 
and the summation is up to a given maximum wavenum- 
ber fc max . Note that the S/N is independent of the bin 
width assumed, as long as the power spectrum does not 
vary rapidly within the bin widths. 

Fig|7] shows the S/N as a function of fc max for the spec- 
tra at z = 0, 1 and 3 in real space (left panel) and in 
redshift space (right), respectively. The results for S/N 
shown here are for a volume of V = l/i~ 3 Gpc 3 (the S/N 
scales with V as S/N oc V 1 / 2 ). For comparison the solid 
curve shows the S/N for the Gaussian error case, which 

scales as S/N oc fcm ax independently of redshift. The 
simulation results show that the non-Gaussian errors de- 
grade the S/N. The degradation becomes more signif- 
icant with increasing fc max and at lower redshifts: for 
the results in real space the S/N is degraded by up to 
a factor of 4 and 2 for z = and 1, respectively, com- 
pared to the Gaussian error case. It should be worth not- 
ing that the S/N becomes nearly constant on fc max ^0.2 
and 0.3/iMpc" 1 , i.e no gain in the S/N even if including 
modes at the larger fc, as has been found in the previous 
works (Rimes & Hamilton 2005, 2006; Hamilton, Rimes 
&: Scoccimarro 2006; Neyrinck, Szapudi & Rimes 2006; 
Neyrinck & Szapudi 2007; Lee & Pen 2008; Angulo et al. 
2008; Smith 2008). 

In the dashed curves we use the perturbation theory 
(PT) in Eq. (O to calculate the covariance, while the 



TABLE 1 
WFMOS Survey Parameters 



Redshift 


Volume 


Number Density 


Bias 




(h- 3 Gpc 3 ) 


(h- 3 Mpc- 3 ) 




0.5 - 1.3 


4.0 


5 X 10- 4 


1.7 


2.3 - 3.3 


1.0 


5 x 10" 4 


3.2 



power spectrum measured from the simulations is used 
for the numerator in the S/N calculation. The PT pro- 
vides a better fit to the data in the linear and weakly 
nonlinear regime, which is coincident within 10% for 
fcmax < 0.16/i/Mpc at z = 0, fc max < 0.23/i/Mpc &tz = l 
and fc max < 0.4/i/Mpc at z = 3, respectively. However, 
at small scale the deviation is so large, since the theory 
predicts the much smaller covariance than the data as 
shown in Figs. [1] and El 

The dotted curves in the left panel of Fig. [7] show the 
halo model results, which fairly well fit the simulation 
results at z = 0. At higher redshifts (z = 1,3), the halo 
model reproduces a saturation in the S/N amplitude on 
the small fc max , but underestimates the impact of the 
non-Gaussian errors, which is due to the underestimation 
in the off-diagonal elements of the covariances in the halo 
model (see Fig. [2]). 

The impact of the non-Gaussian errors on S/N is mit- 
igated in redshift space as in Fig. [3] Also note that 
the S/N in redshift space continues to increase for the 
smaller scales fc max > 0.2/i/Mpc. This is again because 
the nonlinear redshift distortions cause strong suppres- 
sion in the non-Gaussian covariances, making the S/N 
closer to the Gaussian error case. 

We make a more realistic estimate for the S/N tak- 
ing into account the shot noise effect of galaxies that 
are biased tracers of large-scale structure. To do this we 
consider a galaxy survey that resembles the planned sur- 
vey by WFMOS (Wide-Field Fiber-Fed Optical Multi- 
Object Spectrograph), and assume the fiducial survey pa- 
rameters given in Table [T] (see also the WFMOS feasibil- 
ity report 5 ). The target galaxies are supposed to be emis- 
sion line galaxies and Lyman-break galaxies, at z ~ 1 and 
z ~ 3, respectively. To compute the power spectrum and 
covariance of galaxies, we employ a linear bias model for 
simplicity. This assumption is not accurate since the bias 
is generally nonlinear and scale-dependent (Smith et al. 
2007), and hence our results below just give a rough es- 
timate on the S/N. The bias parameters in Table Q] are 
chosen such that the rms density fluctuations of galaxies 
within a sphere of radius 8/i _1 Mpc become a g s = 0.8 
(Glazcbrook et al. 2005). We simply include the effects 
of the linear bias and the shot noise by replacing the 
power spectrum and the covariance in the S/N evalua- 
tion as P — > b 2 P and cov — * 6 4 cov + 2b 2 Pn~ 1 + fig 2 , 
where n g is the mean number density of galaxies. The 
above replacement is done in real space. 

The symbols in Fig. [8] show the simulation results in 
real space. The expected S/N is found to be very signif- 
icant: S/N ~ 400 and 200 for the slices of z = 1 and 3 
for fc max = 0.4/i/Mpc, respectively. This implies that the 
WFMOS-type survey allows a precision of measuring the 

5 Feasibility Study Report : http://www.gemini.edu/files 
/docman/science/aspen/WFMOS-feasibility _report_public.pdf 
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Fig. 7. — The cumulative signal-to-noise ratios are shown as a function of fc max in real space (left) and in redshift space (right), where the 
power spectrum information over 2n/L < k < fcmax is included (L is the box size). Note that the S/N amplitudes are for the simulation 

box volume V = lh~ 3 Gpc 3 . The solid curves in each panel show the S/N for the Gaussian covariance case, which scales as S/N <x fc^ax- 
The simulation results increasingly deviate from the Gaussian error results with increasing fc ma x- In the left panel the dashed curve shows 
the analytical prediction for S/N when the PT trispectrum is used to model the non-Gaussian covariance (the simulation result is used for 
the power spectrum in the S/N calculation). The dotted curves show the results when the halo model is used to model the non-Gaussian 
covariance. 
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kmax (h/Mpc) 

Fig. 8. — The expected S/N for model WFMOS surveys of 
Z ~ 1 and z ~ 3 slices (see Table [1} . Note that the shot noise 
contribution to the covariance is included. The solid curves de- 
note the S/N without the non-Gaussian covariance contribution. 
The bottom panel shows the ratio between the simulation result 
and the solid curve for each redshift slice. The Gaussian error 
assumption overestimates the S/N by 30% (7%) at z = 1(3) for 
fcmax = 0.4/i/Mpc. 



power spectrum amplitudes at a sub-percent level. 6 The 

6 The Fisher information matrix for the 
power spectrum measurement is given as Fu = 
£fc li2 [0P(fcl)/00i] co v - 1 (k 1 ,k2)[dP(k 2 )/de i }, where 0, is a 
set of cosmological parameters of interest. Roughly speaking, 
the unmarzinalized uncertainty in estimating a parameter 0; is 
given as a 2 (6 l ) = [F l ./\- 1 / 2 oc (S/N)- 1 (this exactly holds if 



solid curves are the results obtained assuming the Gaus- 
sian covariances with the shot noise contribution, which 

do not scale as S/N oc fcrnaL on scales where the shot 
noise is relevant in the covariance (n g P ^,1). Compared 
with Fig. one can find that the shot noise causes pos- 
itive and negative effects on S/N: it reduces the overall 
amplitudes of S/N, but mitigates the degradation due to 
the non-Gaussian errors. 

Since the precision of constraining individual cos- 
mological parameters such as dark energy parameters 
roughly scales with the S/N amplitude 6 , Fig. [S] implies 
that the constraining power of the fiducial WFMOS sur- 
vey is degraded by the non-Gaussian errors, compared 
with the Gaussian error case. The impact of the non- 
Gaussian errors on cosmological parameter estimations 
will be presented in a subsequent paper (Takahashi et al. 
in preparation). 

7. EFFECTS OF LONG- WAVELENGTH FLUCTUATIONS 

We have so far employed, as usual, the simulations 
with periodic boundary conditions, where there is no 
clustering power on scales greater than the simulation 
box. Obviously, however, the real universe never obeys 
the periodic boundary condition and does contain the 
density perturbations of scales greater than a surveyed 
volume. In particular, Rimes & Hamilton (2006) pointed 
out a new source of the non-Gaussian errors that in- 
evitably arises when the power spectrum is estimated 
from a finite-size volume, called the beat-coupling (BC) 
effect (2006; also see Hamilton, Rimes & Scoccimarro 
2006; Sefusatti et al. 2006). If the survey region is em- 

the power spectrum amplitude, for the fixed shape, is considered 
as the parameter). Therefore the S/N amplitude gives a rough 
estimate on the precision of parameter estimation provided the 
power spectrum measurement: the greater S/N means the higher 
precision. 
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bedded in a large-scale overdensity or underdensity re- 
gion, then the small scale fluctuations we measure may 
have grown more rapidly or slowly than the ensemble 
average. There are thus non-vanishing correlations of 
the small-scale fluctuations with the unseen large-scale 
fluctuations. This physical correlations may add uncer- 
tainties in measuring the power spectrum on scales of 
interest. 

In this section, therefore, we study how the periodic 
boundary conditions and the density perturbations larger 
than a survey volume (the volume where the Fourier 
transform is performed) affect the power spectrum es- 
timation and the covariance. For this purpose we study 
the following three cases: 

Case 1: We first divide each simulation region of l/i~ 3 Gpc 3 
into eight cubic sub-boxes of equal volume. Each 
sub-box has a volume of (500/i -1 Mpc) 3 and con- 
tains about 128 3 particles. We then randomly se- 
lect only one sub-box and use the particle distri- 
bution to resemble the density perturbation field. 
The density perturbation field outside the sub-box 
is zero-padded within the whole box of l/i~ 3 Gpc 3 . 
The mean mass density is computed from the num- 
ber of the particles within the sub-box. Then we 
perform the FFT of 512 3 grids for the whole box 
to estimate the power spectrum. 

Case 2: Similar to Case 1, but the FFT of 256 3 grids 
is performed only within the sub-box of volume 
(500/i _1 Mpc) 3 that contains the iV-body particles 
(therefore no zero-padded region). 

Case 3: We run new simulations of volume (500/i~ 1 Mpc) 3 
using 128 3 particles and employing the periodic 
boundary condition. Then the power spectrum is 
estimated from the whole box using the FFT of 
256 3 grids. 

Note that the effective mass and spatial resolutions are 
the same in all the cases. We use 400 realizations for 
each case. Cases 1 and 2 do not employ the periodic 
boundary condition and contain the density fluctuations 
larger than the FFT-used volume in structure formation. 
However, these two cases are different in that the funda- 
mental mode of Fourier transform, given as e = 2n/L (L 
is the size of FFT volume), is smaller in Case 1 by factor 2 
than in Case 2. Therefore the density fluctuation field is 
sampled by the finner Fourier modes in Case 1. Also note 
that Case 1 corresponds to a case that the FFT trans- 
form is applied to a survey with complex geometry. Case 
3 has the periodic boundary condition, and is equivalent 
to the procedure we have employed up to the preceding 
section. By comparing these three cases, we will below 
address the effects of the periodic boundary condition, 
the finite Fourier sampling and the beat-coupling effect. 

Fig. [5] shows the real-space power spectra at z = 0, 
estimated according to the procedures described above, 
where the simulated power spectra are divided by the 
non-wiggle linear power spectrum in Eisenstein & Hu 
(1999) for illustrative purpose. The cross, plus and 
dotted symbols are the results for Cases 1, 2 and 
3, respectively. All the results agree well on scales 
k }t O.l/iMpc^ 1 . However, the results for Cases 1 and 
2, which do not impose the periodic boundary condition, 
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Fig. 9. — The real-space power spectra at z = 0, computed 
based on the procedures of Cases 1, 2 and 3 in § [7] Note that 
the power spectra are divided by the non- wiggle linear power spec- 
trum of Eisenstein & Hu (1999) for illustrative purpose. The dot- 
ted symbols show the spectrum estimated from the simulations 
with the periodic boundary condition. The cross and plus symbols 
show the results without the periodic boundary condition: the re- 
sults include contributions from the density perturbations of scales 
greater than the Fourier-transformed volume. The density pertur- 
bation fields for the cross and plus symbols are equivalent, but the 
Fourier-transformed volume for the cross symbol is set to contain 
the zero-padded region (see text for the details). The solid curve 
denotes the linear theory prediction. The dashed curve is the same 
as the solid curve, but convolved with the window function. 

underestimate the power spectrum amplitudes at the lin- 
ear regime k ~ 0.1/iMpc- 1 by up to 10%. This is because 
the non-periodic density fluctuation field is expanded by 
the FFT transform that has periodic basis eigenfunctions 
within the box size (see Sirko 2005 for the similar discus- 
sion). This underestimation can be corrected for if the 
Fourier kernel of the non-periodic field is properly taken 
into account. The dashed curve is the linear power spec- 
trum convolved with the window function which is given 
as W(x) = 1 (= 0) inside (outside) the sub box: 



JV(k) 



l 

V 



r1 z k' ~ 2 

" y P(k') W(k-k') 



(27T) 



The window function in Fourier space is 

sm(fcii/2) 



hL/2 



(13) 



(14) 



with L = 500/i _1 Mpc. The dashed curve is the spherical 
averaged power spectrum, Pwo{k) — J <ifife/(47r) Pvp(k), 
which reproduces the dumping of the power spectrum at 
k < 0.1/i/Mpc. 

FigllOl compares the power spectrum covariances. The 
top two panels show the diagonal parts normalized by 
the linear power spectrum (upper-left) and the spectrum 
measured from simulations (upper-right), which are sim- 
ilar to FigsQ]and[3l respectively. The number of modes 
Nk in the vertical axis is used for the (500Mpc//i) 3 box 
for all the cases. First, comparing the results for Cases 
2 and 3, we find that there are stronger non-Gaussian 
errors at k £ 0.2ft,Mpc -1 for Case 2: the presence of the 
density fluctuations larger than the survey volume, as 
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Fig. 10. — The power spectrum covariances in real space at z = 0, computed from the three numerical settings as in Fig. [9] The diagonal 
components of the covariances as a function of k, normalized by the linear Gaussian covariances (upper-left panel) and by the measured 
Gaussian covariances (upper-right), respectively. The bottom panels show the correlation coefficients r(k\,k2) for the covariances as a 
function of &2, for k\ = 0.035/iMpc — 1 (left) and 0.135/iMpc -1 (right), respectively. In all the panels, the solid curves show the analytical 
predictions for 2 = obtained when the PT is used to model the non-Gaussian covarian ces. The dashed curves show the results when the 
additional non-Gaussian errors due to the long- wavelength fluctuations, modeled by Eg. 1151 . are further included. 




k max (h/Mpc) k max (h/Mpc) k max (h/Mpc) 

Fig. 11. — The signal-to-noise ratios at redshifts z = (left), 1 (center) and 3 (right), computed from the three numerical settings as 
in Fig. [9] The dotted symbols show the equivalent results to Fig. [7] The cross and plus symbols show the results, where the simulations 
without the periodic boundary condition are used. These simulations include contributions from the fluctuations at length scales greater 
than the Fourier volume - the beat-coupling (BC) effect. The solid curves are the result for the Gaussian error case, while the dashed 
curves show the result s of three redshifts obtained when the non-Gaussian covariances are computed from the PT trispectrum plus the BC 
effect modeled by Eq. H15I I. The naive BC model significantly overestimates the non-Gaussian errors seen from the simulations. The dotted 
curves show the result obtained when the halo model trispectrum is further included. 
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in the real universe, increases the non-Gaussian error 
strengths due to the mode-coupling between the large- 
and small-scale density fluctuations. Next, let's compare 
Case 1 with Cases 2 and 3. The results for Case 1 are 
clearly smaller than Cases 2 and 3 even at the linear 
scales such as k 0.1/iMpc^ 1 , but show stronger non- 
Gaussian errors than Case 3 on the large fc's as in Case 

2. The differences between Cases 1 and 2 are caused by 
the presence of the zero-padded regions within the FFT 
volume and the finer Fourier sampling. For Case 1, the 
density perturbations of scales comparable to the FFT 
volume are non-periodic due to a mixture of the zero- 
padded region and the N-body particle distribution in 
the sub-volume. Therefore, the Fourier transform causes 
artificial cross-correlations between the Fourier modes of 
different k's even in the linear regime. Due to the cross- 
correlations, the off-diagonal covariances are amplified as 
shown in the lower panels, while the diagonal covariances 
are relatively suppressed. 

The behaviors of the off-diagonal covariances shown in 
the lower panels are similarly understood. Comparing 
the three cases one can find that the long-wavelength 
fluctuation effect and the zero-padding plus the finer 
Fourier sampling cause stronger cross-correlations be- 
tween the spectra of different /c's over the range of k 
we have considered. 

More important results are given in Figs. [Til showing 
the cumulative S/N values for Cases 1, 2 and 3, com- 
puted properly taking into account the covariances in 
Figs. QUI Each panel shows the results at redshifts of 
z = (left), 1 (center) and 3 (right). First of all, all 
the results agree well with the Gaussian error case on 
the linear scales, k O.l/iMpc -1 . At the larger fc's, the 
results for Cases 1 and 2 show that the presence of long- 
wavelength fluctuations further degrades the S/N ampli- 
tudes by 20% at redshift z — and by 30% at z = 1 and 

3, respectively, compared to the results with the periodic 
boundary condition (Case 3). This implies that, even for 
high redshifts and at the BAO scales, the additional non- 
Gaussian errors due to the long-wavelength fluctuation 
effect need to be included in the analysis for an actual 
survey. Interestingly, the S/N values for Cases 1 and 
2 become to agree well, even though their covariances 
are very different as shown in Fig. 1101 This agreement 
is reasonable, because Cases 1 and 2 contain the similar 
mass density fields of same volume; therefore the amount 
of cosmological information to be extracted is similar. 
These are different only in the FFT procedures. 

A full physical understanding of the complex covari- 
ance behaviors is beyond the scope of this paper. Nev- 
ertheless, it would be interesting to compare the simu- 
lation results with an analytical model. A crude model 
to describe the long-wavelength fluctuation effect on the 
non-Gaussian covariance is proposed in Hamilton et al. 
based on the perturbation theory (also see Takada & Jain 
2008): 

1 /17\ 2 

COV B c(fcl,fe) = yT6 I — J P\in{e)P\in(kl)P\in(k2), 

(15) 

where e = n/L(L = 500/iMpc- 1 for Cases 1 and 2). This 
model ignores the Fourier transform effect of the non- 
periodic density field and rests on a simplified assump- 
tion that the long-wavelength fluctuation effect arises 



from a correlation of the fundamental Fourier mode e 
with the wavenumbers we want to measure. Developing 
a more accurate analytical model of the long-wavelength 
fluctuation effect is now in progress and will be presented 
elsewhere (Kayo et al. in preparation). 

The solid curves in Fig. fTDl show the S/N for the Gaus- 
sian error case, and roughly explains the simulation re- 
sults up to the linear regimes. The dashed curves show 
the PT model predictions including the BC effect mod- 
eled by Eq. (|15p . which are intended to reproduce the 
results for Case 2. The dotted curves show the results 
obtained when the halo model contribution to the covari- 
ance is further included. These analytic models do not 
describe the complex behaviors of the diagonal and off- 
diagonal covariances seen in the simulations. Also the 
simulation results for S/N cannot be explained by the 
analytic models. The analytical model in Eci. tfrS"]) overes- 
timates the beat-coupling effect. This conclusion agrees 
with a recent study of Reid, Spergel & Bode (2008), 
where they showed that the simulation results are fairly 
well explained if Eq. (fT5")) reduced by a factor 3 is added 
to the Gaussian covariance. 

8. DISCUSSION AND CONCLUSION 

Having well-calibrated, accurate covariances of the 
power spectrum is clearly needed in order to obtain 
unbiased, robust cosmological constraints from ongo- 
ing/future BAO experiments. In previous studies, the 
covariance matrix is calculated either by using ana- 
lytic models which cannot be applied to fully nonlinear 
regimes, or by using a limited number of simulation re- 
alizations (Scoccimarro et al. 1999; Meiksin & White 
1999; Rimes & Hamilton 2005 and Neyrinck & Szapudi 
2008). In this paper, we used a very large number (5000) 
of the realizations to study the power spectrum covari- 
ances, allowing us to achieve the convergence at a few % 
level. 

We have carefully studied how the non-Gaussian error 
contributions to the covariance vary with scales and red- 
shifts for the concordance ACDM model. As expected in 
the CDM model, the non-Gaussian errors become more 
significant on smaller length scales and at lower redshifts. 
For redshifts z = 0, 1 and 3, the cumulative signal-to- 
noise (S/N) ratios for measuring the power spectrum 
over 0.01 k 0.4/iMpc _1 are degraded due to the non- 
Gaussian errors by a factor of 1.3, 2.3 and 4, respectively, 
compared to the Gaussian error cases. This degradation 
is slightly mitigated in redshift space because the nonlin- 
ear redshift distortions cause a stronger suppression in 
the covariance amplitudes than in the power spectrum 
amplitudes. 

We also estimated how the density fluctuations of 
scales greater than a survey size cause additional non- 
Gaussian errors via the correlations with the fluctuations 
we want to measure, which inevitably arises for a finite- 
size survey - the so-called beat-coupling effect. This ef- 
fect disappears when estimating the power spectrum co- 
variances from simulations with the periodic boundary 
condition. Thus we rather used the sub-region of the 
original simulation to estimate the new non-Gaussian er- 
rors, and showed that the beat coupling effect can be 
important even in the weakly nonlinear regime and for 
high redshifts: it further suppresses the S/N by 20% at 
z = 0, and 30% at z = 1 and 3, respectively. How- 
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ever, the behaviors of these non-Gaussian errors cannot 
be described by the naive analytic models with and with- 
out the beat-coupling effect. Therefore it will be worth 
exploring a more accurate analytical model of the non- 
Gaussian covariances. Such a model will help us to ob- 
tain physical interpretation and to calibrate the derived 
covariance for arbitrary cosmological models and survey 
parameters (Kayo et al. in preparation). 

We also studied the probability distribution of the 
power spectrum estimators among the 5000 realizations. 
We found that the distribution is nearly Gaussian even 
in the nonlinear regime. More precisely, the mean of the 
power spectrum estimators is not largely biased from the 
ensemble average, and the scatters are well given by the 
diagonal power spectrum covariance at a given wavenum- 
ber. 

A more important question would be how an actual 
galaxy survey is affected by the non-Gaussian errors. For 
this purpose, we made a simplified estimate on the S/N 
expected for a WFMOS-type survey, further taking into 
account the shot noise contamination to the covariance 
due to finite number densities of galaxies. Since the shot 
noise contributes only to the Gaussian errors (in an ideal 
case), including the shot noise not only reduces the to- 
tal S/N amplitude, but also mitigates the influence of 
the non-Gaussian errors. Thus the impact of the non- 
Gaussian errors does vary with survey parameters. Since 
the precision of a given survey for constraining cosmolog- 
ical parameters roughly scales with the S/N amplitude, 
an optimal survey design needs to be realized by taking 
into account the non-Gaussian errors, given the resources 
and observing times for a survey. Furthermore the non- 



Gaussian errors may cause the best-fitting parameters to 
be biased if the model fitting is done improperly assum- 
ing the Gaussian covariances, because the non-Gaussian 
errors are more significant on smaller length scales and 
cause correlated uncertainties between the band pow- 
ers. These issues will be studied in a forthcoming paper 
(Takahashi et al. in preparation). 

Our simulation results of the power spectrum P{k) and 
the covariance matrix cov(fci, k 2 ) are available as numeric 
tables upon request (contact takahasi@a.phys.nagoya- 
u.ac.jp). 
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APPENDIX 

CONVERGENCE OF THE COVARIANCE MATRIX 

It is useful to estimate the necessary number of realizations to achieve a desired accuracy for estimating the power 
spectrum covariance. In this appendix, we examine the numerical convergence of the covariance estimation. Let us 
define the dispersion of the covariance matrix as 

2 , . ((cov(fci,fc 2 ) - (cov(fci,fc 2 ))) 2 ) , A1 s 

CTcovlfcl.tej = -. j- • (Al) 

(cov(«i, k 2 )) z 

Here, (cov(ki,k2)) is the ensemble average over all the realizations, while cov(ki,k2) is the covariance estimated from 
a subset of the realizations whose number is denoted as N r . Fig. [T2l shows the dispersion in Eq. ljAlj) as a function 
of the number of realizations N T . The left panel is for the diagonal elements, and each color symbols correspond to 
fci, 2 = 0.05 (green), 0.2 (blue) and 0.4/i/Mpc (red) with the bin width Ak = 0.01 (circles) and 0.005/i/Mpc (crosses). 
The solid line represents 2/iV r which fits the data very well. Hence, we numerically find the scaling of the dispersion 
given by 

a c 2 ov (fci,fci) ~ (A2) 

Note that the above result is independent of the scale k, the bin width Ak and the simulation box volume. 

The right panel is for the off-diagonal elements for (fci, fc 2 ) = (0.05, 0.2), (0.05, 0.4) and (0.2, 0.4)/i/Mpc, respectively. 
The solid lines represent 10/Ar and 100/A r . Similar to the diagonal parts, we obtain 

c 2 ,(fci , fc 2 ) oc — — — -77 cx . (A3) 

covV ' ' N r (Ak) 2 N r AN kl AN k2 y ' 

for k\ ^ fc 2 . Here AN ki (i — 1,2) are the numbers of modes available for the bins fej with the bin width, and the 
proportionality factor would depend on the scale fcj. The analytical derivation for Eqs. (|A2p and (|A3[) will be presented 
in Kayo et al. (in preparation). 

PROBABILITY DISTRIBUTION OF THE POWER SPECTRUM ESTIMATOR IN GAUSSIAN LIMIT 

In the linear regime, the real and imaginary parts of the density fluctuation (Re[<5k] an d Im[<5k]) follow the Gaussian 
distribution with mean and dispersion P(k)/2. The power spectrum estimator for a given realization P is the 
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diagonal off diagonal 




number of realizations N r number of realizations N r 

Fig. 12. — The dispersions among the power spectrum covariances each of which i s est imated from the N T realizations (a subset of the 
while 5000 realizations), as a function of N T . The dispersion is estimated using Eq. jAlt . The left (right) panel show the results for the 
diagonal (off-diagonal) parts for varying the wavenumber bins and the bin widths. The color symbols are the simulation results, while the 
solid curves denote the approximate fittings (see text for the details). The plots explicitly show that the power spectrum covariances are 
estimated at a sub-percent level accuracy by using our whole 5000 realizations. 



summation of the squared Gaussian fluctuations, P(k) = (1/^Vfe) X)k(^- e [^k] 2 + Iui[<5k] 2 )- Then the distribution of P 
obeys the chi-square distribution (e.g. Abramowitz & Stegun 1970): 



F(P(k);N k /2) = * . . ( ^t^M e -Pikyp(k) 



N k /2 



T(N k /2) \ 2 P{k) I p(jfc) 



(Bl) 



for P > and F = for P < 0. Its mean and dispersion are P and P 2 / (Nk/2), respectively. The skewness and 
kurtosis are given in Eq. ljlip . In the limit of N k — > oo it reduces to the Gaussian distribution. The factor two in the 
number of modes N k /2 arises because 5^ and are not independent. 
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